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Abstract 

The recently developed Meron-Cluster algorithm completely solves 
the exponentially difficult sign problem for a number of models previ- 
ously inaccessible to numerical simulation. We use this algorithm in a 
high-precision study of a model of N = 1 flavor of staggered fermions in 
(2+l)-dimensions with a four-fermion interaction. This model cannot be 
explored using standard algorithms. We find that the Z(2) chiral sym- 
metry of this model is spontaneously broken at low temperatures and 
that the finite-temperature chiral phase transition is in the universality 
class of the 2-d Ising model, as expected. 
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1 Introduction 



There are a number of models of interest which suffer from a very severe sign prob- 
lem. This includes QCD and other field theories with a non-zero chemical po- 
tential or a non-zero vacuum angle or odd numbers of fermion flavors, frustrated 
quantum spin systems, like the quantum antiferromagnet in an external magnetic 
field, and models for strongly-correlated electrons, like the Hubbard model for high- 
temperature superconductivity. These models have a Boltzmann weight which can 
be negative or even complex and so cannot be interpreted as a probability. This 
difficulty can be overcome in numerical simulations by including the sign or phase 
of the Boltzmann weight with observables. Unfortunately, this leads to large cancel- 
lations and gives exponentially small observables. This requires exponentially large 
statistics, which makes it in practice impossible to simulate these models numeri- 
cally. 

Recently, a new technique has been developed, called Meron-Cluster algorithms 
fl|], which completely solves the sign problem for some of these models [@, f|. It 
identifies the origin of the sign problem with properties of the clusters, which enables 
it to be eliminated. Cluster algorithms in general are extremely efficient at exploring 
configuration spaces and very often do not suffer from critical slowing down as a 
phase transition is approached, unlike many other algorithms. For example, in this 
and previous papers, we can work directly in the chiral limit with massless fermions. 
Combined with the ability to construct improved estimators, we can perform a high 
precision study of these models with only modest statistics. 

In this paper, we explore a model of N — 1 flavor of staggered fermions in 
(2+l)-dimensions with a four-fermion interaction. This model has a very severe 
sign problem and cannot be simulated with standard techniques. We build a Meron- 
Cluster algorithm, which we use to perform a high-precision study. We find that the 
Z(2) chiral symmetry of this model is spontaneously broken at low temperatures 
and, using finite-size scaling analysis, we verify that the finite-temperature chiral 
phase transition is in the universality class of the 2-d Ising model. A recent study of 
the same model with Z(2) chiral symmetry in (3+l)-dimensions has shown, using 
a Meron-Cluster algorithm, that a finite temperature chiral phase transition occurs 
which has the universal behavior of the 3-d Ising model |J. The work presented in 
this paper concerns a different universality class and also constructs more observables 
than were previously considered. 

The identification of the finite temperature critical behavior is not entirely straight- 
forward 0. A model of N fermion flavors with a four-fermion interaction shows 
mean-field behavior in the N = oo limit. On the other hand, at finite N one finds 
the non-trivial critical behavior that one expects based on dimensional reduction 
and standard universality arguments. For example, in |J it has been verified that 
the chiral phase transition in a (2 + l)-d four-fermion interaction model with N = 4 
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flavors and Z(2) chiral symmetry is in the universality class of the 2-d Ising model. 
Due to the fermion sign problem, standard fermion simulation methods often do 
not work in models with too small a number of flavors. The work presented in this 
paper shows that the same universal behavior holds for N — 1 flavor. 

The standard technique to deal with fermions in Monte-Carlo simulations is to 
integrate them out, resulting in a non-local bosonic theory. The Meron-Cluster al- 
gorithm does not integrate out the fermions, instead we describe them with a local 
theory using a Fock space basis of occupation number. The fermion sign arises as a 
non-local property due to the permutation of fermion world lines. Using probabilis- 
tic rules, we connect neighboring lattice sites, producing closed loops, which are the 
clusters. A cluster is flipped by making all of its occupied sites empty and its empty 
ones occupied. Such a cluster flip can change the fermion sign by changing the 
permutation of fermion world lines. A cluster whose flip changes the sign we call a 
meron. We can tell if a cluster is a meron simply from its structure. A typical config- 
uration contains many merons, yet the observables of interest are only non-zero for 
configurations with very few or no merons. The signals from standard Monte Carlo 
algorithms are so exponentially small because the Markov chain explores a vast con- 
figuration space, yet only an exponentially small sub-space makes any contribution 
to measurables. By restricting ourselves to only explore the relevant sub-space, we 
completely solve the sign problem. 

This paper is organized as follows. In section 2, we present the fermionic model 
which we have studied, calculate its partition function using the Hamiltonian for- 
mulation and find that there is a sign problem. In section 3, we describe briefly 
the Meron-Cluster algorithm which we have used to perform numerical simulations 
of this model. We present the results of the simulations in section 4 and give our 
conclusions in section 5. 



2 The Staggered Fermion Model 

We consider staggered fermions in the Hamiltonian formulation on a 2-dimensional 
spatial lattice of extent L, which is even. The Hamiltonian operator is 

x,i x 

h x , = ifeiW^ + ^^ + GC^.-i)^^-!), (2.1) 

where r] Xt i = 1 and T] Xj 2 = {— are the standard Kawamoto-Smit phases for 
staggered fermions and G is a constant. The fermionic operators satisfy the usual 
anticommutation relations {^f x ,^ y } = {^J, = 0, {^^,^ y } = 5 xy . The same 
model in (3+l)-dimensions was explored in [f|. We refer the reader to this paper, 
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where various features of the model and the Meron-Cluster algorithm are discussed 
in more detail than we give here. 

In the Hamiltonian formulation of the theory, fermion doubling on the lattice 
occurs only in the spatial dimensions. Using staggered fermions, the Dirac compo- 
nents of a spinor are distributed spatially, reducing the number of fermion flavors 
by a factor of four. Thus this (2+l)-dimensional model contains N — 1 fermion 
flavor. The model has a global U(l) symmetry corresponding to conserved particle 
number, as the total particle number operator commutes with the Hamiltonian 

N = J2*5*x, [H,N]=0. (2.2) 

X 

Furthermore the Hamiltonian has, for m = 0, a discrete Z(2) symmetry correspond- 
ing to shifts by one lattice spacing. However the mass term breaks that symmetry 
explicitly. For a single flavor of massless fermions, the symmetry of the lattice model 
is U(l) ® Z(2) and we refer to the discrete symmetry as chiral symmetry. In the 
continuum, a single massless fermion flavor has a U(l) axial symmetry (there is 
no gauge interaction, so this symmetry is not anomalously broken). The discrete 
Z(2) symmetry is the lattice remnant of this continuous symmetry. From now on, 
we set m = and explore the behavior of the chiral symmetry of this model. The 
symmetries of staggered fermions are discussed in detail in Ref. 0. If the Z(2) 
chiral symmetry is spontaneously broken at some finite temperature, from univer- 
sality we expect this to be a second-order phase transition. As the critical point is 
approached, the correlation length £ diverges and the system becomes insensitive 
to the time extent. Due to dimensional reduction, we expect a finite-temperature 
chiral phase transition in this model to belong to the 2-d Ising universality class. 

The partition function of the model is 

Z = Tr[exp(-/3F)] = lim Tr[exp(-eiJ)] M (2.3) 

M— >oo 

= lim Tr[exp(— eHi) exp(— eH 2 ) exp(— eif 3 ) exp(— eH 4 )] M , 

where we use the Suzuki- Trotter decomposition to divide the Euclidean time extent 
(3 into 4M time slices, the lattice spacing in the time direction being e = j3/M. 
The Hamiltonian operator is decomposed into four parts H — H\ + Hi + H3 + H4. 
All of the terms that contribute to a particular Hi commute with one another, as 
each term is an interaction between nearest-neighbors and each lattice site appears 
in only one such nearest-neighbor pair. However, the Hi do not commute with one 
another. We note that it is not actually necessary to discretize the time direction, 
as it is possible to work directly in the Euclidean time continuum 0]. 

We can equivalently describe this model with bosonic operators, using a trans- 
formation by Jordan and Wigner [jSJ. We order the lattice sites on each time slice 
arbitrarily into a chain, which can be done in any number of spatial dimensions. 
For example, a possible ordering of points in two spatial dimensions is by an index 
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I — x\ + (x2 — l)L. The fermionic operators are now represented by a chain of Pauli 
matrices 



= o*<4...of_ 1 v+, tt x = 0*0*. = + 1) (2-4) 
a ± = ^(a 1 ±ia 2 ), [<r>4] = 2tf, m e*M, 

where the spatial position x is denoted by the index / and the Pauli matrices satisfy 
the usual commutation relations. To calculate the partition function of the theory, 
we use the Fock space basis of occupation number n x — 0, 1 i.e. the eigenstates 
of a 3 . The occupied and empty states are respectively |1) and |0), which satisfy 
( x 3 |l) = |l)anda 3 |0) = -|0). 

The time evolution operator exp(— eifj) acts on a time slice of occupation num- 
ber states, producing the next time slice. This is decomposed into the product of 
operators exp(— eh Xji ) acting on nearest-neighbor occupation states. The transfer 
matrix is 



exp(-e/i x>i ) = exp(^) 



/ exp(-f ) \ 

cosh I E sinh | 



S sinh I cosh | 



\ exp(-f ) / 



(2.5) 
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the basis being 1 00) , |01), 1 10) and |11), where e.g. 1 1) represents state |0) at x and 
|1) at x + i. If these nearest-neighbors are labelled / and m, the off-diagonal transfer 
matrix elements have a factor E = ^r,j^f+i°f+2"-°"m-i- Note that this operator is 
diagonal in the occupation number basis. 

The partition function of the theory is given as a path integral 

Z f = J2 Sign[n] exp(-S[n]), (2.6) 

n 

where we sum over all possible configurations of occupation numbers n(x, t) — 0, 1 
on a (2 + l)-d space-time lattice of points (x,t). The Boltzmann factor exp(— S[n\) 
for a configuration is the product of the Boltzmann factors for each space-time 
plaquette exp(— s[n(x, t),n(x + i, t),n(x, t + l),n(x + i,t+ 1)]), which are 

eG 

exp(-s[0, 0, 0, 0]) = exp(-s[l, 1,1,1]) = exp(-— ), 

exp(-s[0, 1, 0, 1]) = exp(-s[l, 0, 1,0]) = cosh |, 

exp(-s[0, 1, 1, 0]) = exp(-s[l, 0, 0, 1]) = sinh ^. (2.7) 

All other plaquettes are illegal and have Boltzmann weight zero, as they repre- 
sent non-conservation of fermion number. Any configuration which contains illegal 
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plaquettes has itself Boltzmann weight zero and makes no contribution to the par- 
tition function. We are only interested in legal configurations, which have to satisfy 
several constraints. Note that here we have dropped the overall factor exp(eG/4) 
that appeared in eq. (|2.5|) . The sign of a configuration, Signfn], is also a product of 
space-time plaquette contributions sign[n(x, t), n(x +i,t), n(x, t + 1), n(x + i,t+ 1)] 
with 

sign[0, 0, 0, 0] = sign[0, 1, 0, 1] = sign[l, 0, 1, 0] = sign[l, 1, 1, 1] = 1, 
sign[0,l,l,0] = sign[l,0,0,l] = S. (2.8) 

The occupied lattice sites define world-lines of fermions, which close due to the 
periodicity of the Euclidean time direction. The world-lines are free to permute 
during their time evolution as the fermions interchange position and each configura- 
tion has a well-defined permutation of fermions. The Pauli exclusion principle tells 
us that the sign of a configuration is the permutation sign of the fermions, hence 
Sign[n] = ±1. This non-local effect is contained in the factors £ of each space-time 
plaquette. 

The expectation value of a fermionic observable A[n] is given by 

(A), = ^£A[n]Sign[n]exp(-S[n]) = (2.9) 

(Sign) = — ^Sign[n]exp(-S'[n]), 
z & n 

where (...) means a measurement made in the bosonic ensemble, whose partition 
function is Z b = J2n ex P(~ S[n}). To measure one fermionic observable requires two 
bosonic measurements. The quantities of physical interest which we measure are 
the chiral condensate the chiral susceptibility x an d a Binder cumulant U of 
the chiral condensate, respectively 

^*[n] = 7£(-ir^WM)-5), 

4 x,t Z 

x = JL/myf) u=i — . (2.io) 

X (3V KK ' ;/ ' 3[((^) 2 )/] 2 1 J 



3 The Meron- Cluster Algorithm 

We now describe briefly the Meron-Cluster algorithm which we used to sample the 
bosonic ensemble corresponding to the fermionic model without the sign factor. We 
set G = 1, for which the bosonic model is the isotropic antiferromagnetic quantum 
Heisenberg model, whose Hamiltonian is H = Y)_ ASlS 1 , ~. + S%S 2 , ~. + S?.S 3 , 0, where 

O t—IX,l\ X x+l x X+l x X + l' 1 

S l x = \a\ is a spin 1/2 operator at the lattice site x, labelled by I in the Jordan- 
Wigner chain. There already exist extremely efficient cluster algorithms to simulate 
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weight 



exp (— 7 



cosh (f ) 



sinh (§) 



configuration 



break-ups 



I I 

1 



I I 



V 



1 — p 



1 



Table 1: Cluster break-ups of various plaquette configurations together with their 
probabilities, where p = 2/[l + exp(e/2)]. The dots represent occupied sites and the 
fat lines are the cluster connections. 



bosonic quantum spin systems [11], [12]| , and the first cluster algorithm for lattice 



fermions was constructed in JT3|. These algorithms can be implemented directly 
in the time continuum J7j, i.e. the Suzuki- Trotter time discretization is not even 
necessary. In this study, we discretize the time direction. 

We use the same algorithm that was used in 0. Each configuration is de- 
composed into a set of clusters, which consist of connected lattice sites. A new 
configuration is generated by flipping the clusters. When a cluster is flipped, all 
lattice sites contained in that cluster change occupation number from n(x, t) to 
1 — n(x, t), i.e. the occupied sites become empty and the empty ones occupied. To 
build the clusters, a probabilistic choice is made in each space-time interaction pla- 
quette [n(x, t),n(x + i, t),n(x, t + 1), n(x + i,t + 1)] as to which neighboring lattice 
sites are connected to one another. A cluster is a sequence of connected sites. In 
this algorithm, the clusters are closed loops. The probabilistic choices (called cluster 
break-ups) which build the clusters are designed to obey detailed balance and we 
only allow break-ups which generate legal plaquettes under cluster flips. The cluster 
rules are illustrated in Table 1. For plaquette configurations [0, 0, 0, 0] and [1, 1, 1, 1], 
i.e. entirely empty or entirely occupied, we always connect sites with their time-like 
neighbors. For configurations [1,0,0,1] and [0,1,1,0] where a fermion hops to a 
neighboring site, we always connect sites with their space-like neighbors. For con- 
figurations [1, 0, 1, 0] and [0, 1, 0, 1], i.e. a static fermion next to an empty site, we 
connect the sites with their time-like neighbors with probability p = 2/[l + exp(e/2)] 
and with their space-like neighbors with probability 1 — p. This algorithm was also 
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Sign[n] = 1 



Sign[n] = — 1 




Figure 1: Two configurations of fermion occupation numbers in (1 + 1) dimensions. 
The dots represent occupied sites and the shaded plaquettes carry the interaction. 
With periodic spatial boundary conditions, the two fermions interchange positions 
in the second configuration, giving it Sign[n] = — 1. Flipping the meron-cluster 
(represented by the fat line) changes one configuration into the other, changing the 
fermion sign. 



used in JT^]. It is extremely efficient, has almost no detectable autocorrelations and 
its dynamical exponent for critical slowing down is compatible with zero. 

Each cluster has two orientations, with lattice site occupancies n(x, t) and 1 — 
n(x, t). When a cluster is flipped, the new configuration which is generated may have 
a different sign from the previous one, depending on whether or not the permutation 
of fermion world-lines is changed. A cluster whose flip changes Signfra] we call a 
meron, those which leave Sign[n] unchanged we call non-merons. Flipping a meron 
changes the topology of the fermion world-lines. The term meron has been used 



before to denote half-instantons such as in the 2-d 0(3) model at non-zero 
vacuum angle 9 [ I6[ . The number of merons in a configuration is always even, as 
flipping all clusters leaves the sign unchanged. An example of a meron-cluster is 
given in Figure 1. When the meron-cluster is flipped the first configuration with 
Sign[n] = 1 turns into the second configuration with Sign[n] = —1. For cluster 
algorithms more general than the one described here, it is not always possible to 
identify certain clusters as merons fl7| . 

The meron concept alone gives us an exponential gain in statistics. Starting from 
a configuration containing Nc clusters, we consider the ensemble of 2 N ° configura- 
tions where we allow all possible cluster orientations. If a configuration contains no 
merons, all configurations in the ensemble have Sign[n] = 1. However, if it contains 
merons, half the ensemble has Signfn] = 1 and the other half Sign[n] = — 1, which ex- 
actly cancel, giving a contribution 0. The improved estimator gives (Sign) = (5n,o), 
i.e. the probability that a configuration contains N = merons, which is an expo- 
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nential improvement on standard algorithms, which measure a statistical average of 
±1. As explained in H, this solves half of the sign problem. 



We also construct improved estimators for observables. The chiral susceptibility 



is 



X 



— ((W) 2 )f 



1 ((^) 2 Sign) 
W (Sign) 



(3.1) 



The total chiral condensate for a given configuration, ^^[n] = J2c ^^c, is a sum 
of cluster contributions. Averaging \ over the ensemble of 2 Nc configurations gives 



A' 



(E 



c 



|^^c| 2 5tv,o + 2|^^ Ci 



\^C 2 \$N,2) 



(3.2) 



This only gets contributions from configurations with = or = 2 merons (C\ 
and C*2 are the two merons). The vast majority of configurations contain many 
merons, but they make no contribution to observables. The zero- and two-meron 
sectors of configuration space are exponentially small, but they contain all of the 
contributions to x- Restricting ourselves to only explore this sub-space, we expo- 
nentially enhance both the numerator and denominator of eq. ( |3.2j ), leaving the ratio 
invariant. This solves the remaining half of the sign problem. 

For the Binder cumulant U, we need to measure ((\E'\I/) 4 )/ and hence 

(Sign(^) 4 ) = (Sign ]T ^ c m c m Ck m Cl }. (3.3) 

A cluster's condensate contribution ty^c changes sign when the cluster is flipped. 
When a meron-cluster is flipped, Sign is changed. The non-zero terms in (Sign^^) 4 ) 
do not change sign if any cluster in the configuration is flipped. These non-zero 
terms must contain odd powers of ty^c for all merons C in the configuration and 
even powers of ^^c' f° r ah non-merons C . The average over the ensemble of 2 Nc 
configurations is 



(Sign(^) 4 ) 2 iv c = 5 Nfi 



l^cl 4 + 6 J2 \^c\ 2 \^c\ 2 



(3.4) 



+0~N,2 

+5na 



4|^^ Cl | 3 |^^c- 2 | + 4|^ C2 | 3 |^ Cl | + 12^ \^^ c \ 2 \^^ Cl \\^^ C2 \ 



c 



24| 1 1 W C2 1 1 **c 3 1 1 ^^c 4 1 



where N is the number of merons in the configurations, Ci, 62,63 and C4 are the 
merons and all sums in eq. ( |3.4|) are over non-meron clusters. This average only 
gets contributions from the zero-, two- and four-meron sectors and so we need only 
explore this sub-space. We average this quantity over the complete bosonic ensemble 
to measure (Sign(\I>\I/) 4 ) and hence U. 
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Consider the case of measuring \- We expect that p(0)/p(2) oc (\C\/V/3) 2 , where 
p(0) and p(2) are the probabilities that a configuration has zero or two merons and 
|C| is the average cluster size. In large volumes, the majority of configurations 
has two merons, contributing to (Sign). For even greater accuracy, we reweight 
the meron-sectors with trial probabilities pt{0) and Pt{2), so that they appear with 
roughly equal frequency. This gives 

= (Eg \VVc\ 2 8n,o Pt{0) + 2\^ Cl \\^c 2 \S N ,2 p t {2)) ( v 

X /3V(6 N>0 Pt (0)) ' 1 ' ' 

The reweighting probabilities can be adjusted to minimize the statistical error. This 
technique was previously used in [Tj|. To measure the Binder cumulant U, we use 
reweighting probabilities p t (0),p t (2) and Pt(4). 



4 Numerical Results 



We have performed simulations of the staggered fermion model on lattices with 
antiperiodic spatial boundary conditions from L = 4 up to L = 30 at inverse tem- 
peratures in the range (3 G [1.0,3.0], which includes the critical temperature where 
the chiral symmetry is spontaneously broken. We have made separate runs with 
either a fixed number of time slices (typically M = 10, i.e. 40 time slices) or with 
fixed lattice spacing in the time direction (e = 0.1). In each simulation, we have 
made at least 1000 thermalization sweeps followed by 10000 measurements, with 
these numbers increased by a factor of 10 for L < 10. In one sweep of the lattice, 
a new cluster connection is proposed on each interaction plaquette and each cluster 
is flipped with probability 1/2. To find the optimal reweighting probabilities pt(N) 
which minimize the statistical error, we first make a sample run without reweighting, 
only exploring the relevant meron-sectors. The observed relative weights are then 
used in production runs, where the sectors appear with equal probability. The major 
part of the sign problem is removed by the improved estimators, but reweighting is 
necessary for accurate measurements in large volumes. 

A sample of the data measured is given in Table 2. The Table contains (Sign) 
and the susceptibility \ measured over all meron-sectors, and the reweighted (Sign) r 
and Xr measured over the zero- and two-meron sectors only with the reweighting 
factor p t (0)/p t (2). All of these data are produced with 1000 thermalization sweeps 
and 10000 measurements. As all of the contributions to \ come from the zero- 
and two-meron sectors, % and Xr should be identical. Note that (Sign) r , the frac- 
tion of zero-meron configurations generated by sampling the zero- and two-meron 
sectors only, is typically a lot bigger than (Sign), the fraction of zero-meron con- 
figurations generated over all meron sectors. In small space-time volumes, x can 
be accurately measured even when sampling all meron sectors. However, in large 
space-time volumes, (Sign) is too small to be measured and we can only determine 
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L 


P 


(Sign) 


X 


Pt(0)M(2) 


(Sign) r 


Xr 


8 


1.0 


0.804(5) 


0.829(3) 


0.5/0.5 


0.820(7) 


0.826(6) 


8 


1.5 


0.465(9) 


2.84(3) 


0.5/0.5 


0.52(1) 


2.84(4) 


8 


2.0 


0.214(6) 


9.2(2) 


0.3/0.7 


0.474(9) 


9.0(1) 


8 


2.4 


0.140(4) 


16.6(3) 


0.2/0.8 


0.501(9) 


16.4(3) 


10 


2.4 


0.057(3) 


24.8(6) 


0.2/0.8 


0.369(8) 


24.2(5) 


12 


2.4 


0.0203(8) 


33(1) 


0.1/0.9 


0.443(7) 


34.0(7) 


14 


2.4 


0.0052(6) 


41(4) 


0.1/0.9 


0.338(8) 


44(1) 


16 


2.4 


0.0005(2) 


80(40) 


0.075/0.925 


0.314(4) 


57(1) 


20 


2.4 






0.03/0.97 


0.355(9) 


82(3) 


24 


2.4 






0.01/0.99 


0.46(1) 


120(5) 


28 


2.4 






0.01/0.99 


0.329(9) 


156(8) 



Table 2: Numerical results for the non-reweighted (Sign) and susceptibility x mea- 
sured over all meron sectors, and the reweighted {Sign) r and Xr measured over the 
zero- and two-meron sectors only with a reweighting factor p t (0)/p t (2). For the 
larger volumes, (Sign) and x cannot be measured. 
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Figure 2: The meron number probability distribution for various spatial sizes L 
8,12,20 and 28 at (3 = 2.4. 
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Figure 3: The chiral susceptibility x (is a function of the inverse temperature f3 for 
various spatial sizes L = 8, 10, 12, 14 and 16 on lattices with 40 time slices. The 
chiral symmetry is intact at small f3 and spontaneously broken at large (5. 



the susceptibility by restricting ourselves to the zero- and two-meron sectors. The 
staggered fermion model suffers from a very severe sign problem which is solved by 
the Meron-Cluster algorithm. 

Figure 2 shows the meron number probability distribution in an algorithm that 
samples all meron sectors without reweighting. For small volumes the zero-meron 
sector and hence (Sign) are relatively large, while multi-meron configurations are 
rare. On the other hand, in larger volumes the vast majority of configurations 
has a large number of merons and hence (Sign) is exponentially small. For ex- 
ample, an extrapolation from smaller volumes gives a rough estimate for the non- 
reweighted (Sign) 10~ 9 on the L = 28 lattice at f3 — 2.4., while the reweighted 
(Sign) r = 0.329(9). Even if the configurations are entirely uncorrected, to achieve 
a similar accuracy without the meron-cluster algorithm one would have to increase 
the statistics by a factor 10 18 , which is obviously impossible. In fact, at present 
there is no other method that can be used to simulate this model. 

Figure 3 shows the chiral susceptibility x as a function of f3 for various spatial 
sizes L. At high temperatures (small j3) x is almost independent of the volume, 
indicating that chiral symmetry is intact. On the other hand, at low tempera- 
tures (large (3) x increases with the volume, which implies that chiral symmetry is 
spontaneously broken. To study the critical behavior in detail, we have performed 
a finite-size scaling analysis for x focusing on the range (3 G [2.2, 2.6] around the 
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Figure 4: T7ie chiral susceptibility plotted against L for various values of (3 computed 
with e = 0.1. The fit is to the finite-size scaling formula (U-A) with a(x) expanded to 
first order and b(y) expanded to third order. The exponents are set to the 2-d Ising 
model values. All curves are obtained from one fit. The x 2 V er degree of freedom is 
0. 84 indicating a good agreement of our data with the finite-size scaling ansatz and 
2-d Ising model critical exponents. 



critical point. Since a Z(2) chiral symmetry is spontaneously broken at finite tem- 
perature in this (2 + l)-d model, one expects to find the critical behavior of the 2-d 
Ising model. The corresponding finite-size scaling formula valid close to (3 C is 



X (L,(3)=a(x)+b(y)L^, 



c; 



a{x) = ao + a,\x + a 2 x + x = j3 — (3, 
b(y) = b + b lV + b 2 y 2 + ...,y = (J3- (i c )L Xlv - (4-1) 

For the 2-d Ising model the critical exponents are given by v = 1.0 and 7/1/ = 1.75. 
Assuming these values for the exponents, we obtain (3 C = 2.43(1) for fixed e = 0.1 
from the finite-size scaling fit, with a chi squared per degree of freedom of 0.84. The 
fit of the data is plotted in Figure 4. The value of (3 C is slightly dependent on e. 

In the finite-size scaling equation ( |4.1|) , for large enough L one can neglect the 
term a(x). Then x/L 7/V is a function of y = ({3—f3 c )L 1 / l/ alone, i.e. the susceptibility 
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Figure 5: Finite-size scaling behavior of the chiral susceptibility \. The data for 
various spatial sizes L and inverse temperatures [3 fall on one universal curve. 



data in various volumes at various j3 can be described by one universal function. We 
have varied the value f3 c to find if all the data can be collapsed onto one universal 
curve. In Figure 5, we plot the universal curve obtained by taking j3 c = 2.43. The 
excellent agreement over a large range of spatial volumes L and inverse temperatures 
(5 is an indication of the quality of the finite-size scaling fit. 

We also measure Ul, the Binder cumulant in volumes of extent L. In Figure 6, 
we plot the expected behavior of Ul as L increases for different temperatures. For 
T > T c , the chiral symmetry is intact and U L flows into the T = oo fixed point 
U = 0. For T < T c , the chiral symmetry is spontaneously broken and Ul flows 
into the T = fixed point U = 2/3. If the universality class has a non-trivial fixed 
point U = £/*, then Ul flows into this value at T = T c . By measuring Ul in various 
volumes at many different temperatures, we determine this flow numerically. We 
have measured the Binder cumulant values in volumes up to L = 30 and we plot 
some of these values as a function of 1/L in Figure 7. These measurements are made 
with the number of time slices fixed at 40. Each curve in the figure represents some 
fixed temperature. In Figure 7, for small j3 (i.e. high temperatures), Ul clearly flows 
into the infinite temperature fixed point U — 0, while for (3 large (low temperatures), 
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Figure 6: The expected flow of Ul as a function ofl/L. On each curve the temper- 
ature is constant. 
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1/L 

Figure 7: Measured values of Ul plotted versus 1/L for various f3 in volumes with 
40 time slices. Near f3 = 2.35, the values appear to flow into the non-trivial fixed 
point U* = 0.60(1). 
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Ul flows into the zero temperature fixed point U = 2/3. For (3 close to (3 C , we have 
to go to larger volumes to see this behavior. Near (3 = 2.35, the cumulant values 
appear to flow into a non-trivial fixed point U*. Examining this region closely, we 
estimate the critical inverse temperature as (3 C = 2.36(2) and the fixed point value 
£/* = 0.60(1). The finite-size scaling fit of \ measured at this e ~ 0.24 gives the same 
value of f3 c . Note that this deviates slightly from the critical temperature measured 
at e = 0.1. The universal fixed point value for the 2-d Ising model is estimated as 
£/* ~ 0.58 ||19|| . This is further evidence that the chiral phase transition belongs to 
the 2-d Ising universality class. 



5 Conclusions 



The Meron-Cluster algorithm has recently been developed to allow numerical sim- 
ulations in models which suffer from a very severe sign problem. In this paper, 
we have applied this technique to investigate a model of staggered fermions. Unlike 
standard methods, which integrate out the fermions, resulting in a non-local bosonic 
action, we use a Fock space of occupation number to describe the fermions. We have 
a local bosonic action, with an additional non-local sign factor which contains the 
Fermi statistics. Due to the Pauli exclusion principle, configurations which have an 
odd permutation of fermion world lines have a negative sign. This sign leads to very 
large cancellations in observables and usually makes it impossible to make accurate 
measurements in numerical simulations. The Meron-Cluster algorithm decomposes 
every configuration into closed loops of connected sites, each loop being a cluster. 
Loops which change the fermion sign when flipped are identified as meron-clusters. 
A meron-cluster identifies a pair of configurations with equal weight and opposite 
sign. This results in an exact cancellation of two contributions ±1 to the path 
integral, such that only configurations without merons contribute to the partition 
function. Observables only receive contributions from configurations which contain 
very few or no merons, whereas the vast majority of configurations contain many 
merons. By only exploring the sectors of configuration space with the relevant num- 
bers of merons, one makes an exponential gain in statistics. Combined with efficient 
re-weighting of the remaining meron sectors, this completely solves the sign prob- 
lem. Cluster algorithms are extremely efficient at exploring configuration spaces 
and generating uncorrelated configurations and generally do not suffer from critical 
slowing down. Even in models without a sign problem, the Meron-Cluster algorithm 
is more efficient than standard fermion simulation methods. 

In this paper, we examined a model of N — 1 flavor of staggered fermions 
in (2 + l)-dimensions, which has a Z(2) chiral symmetry. The model has a very 
severe sign problem and cannot be solved by standard fermion simulation algorithms. 
Using a Meron-Cluster algorithm, we were able to make high-precision measurements 
of the chiral susceptibility and Binder cumulant even in very large volumes and 
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low temperatures. In order to perform an accurate and reliable finite-size scaling 
analysis, it was necessary to go to volumes so large, where the sign problem is so 
severe, that a standard algorithm would require statistics on the order of 10 18 to 
attain a similar accuracy. We were able to verify that the model undergoes a finite- 
temperature chiral phase transition, which belongs to the universality class of the 
2— d Ising model. This is the behavior expected from dimensional reduction and 
universality. The same universal behavior was observed in the N = 4 flavor case 
IJ. However, the standard fermion algorithm that was used in that study does not 
work for N < 4 due to the fermion sign problem. 

It is quite natural to use cluster algorithms in models of discrete variables. A 
future possible application of the Meron-Cluster algorithm is in exploring quantum 
link models [2(| which are used in the D-theory formulation of QCD [21, 23|. In 



D-theory, a model of discrete quantum variables undergoes dimensional reduction, 
resulting in an effective theory of continuous classical variables. In quantum link 
QCD the quarks arise as domain wall fermions. The application of meron-cluster 
algorithms to domain wall fermions is in progress. Also there are many applications 
to sign problems in condensed matter physics. Investigations of antiferromagnets in 
a magnetic field and of systems in the Hubbard model family are given in Ref.0. 

At present, the Meron-Cluster algorithm is the only method that allows us to 
solve the fermion sign problem. A severe sign problem arises in lattice QCD calcu- 
lations at non-zero baryon number due to a complex action. It is therefore natural 
to ask if our algorithm can be applied to this case. At non-zero chemical potential 
the 2-d 0(3) model, which is a toy model for QCD, also suffers from a sign problem 
due to a complex action. When applied to the D-theory formulation of this model, 
the Meron-Cluster algorithm solves the sign problem completely ||. It is an open 
question if such progress can be made in investigations of QCD. 
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